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Abstract 


The  objective  of  this  program  was  to  develop  an  accurate  and  fast  simulation 
methodology  for  supercritical  injection  and  mixing  using  molecular  dynamics  (MD).  The  use  of 
MD  allows  the  supercritical  mixing  to  occur  ‘naturally’  without  the  need  for  tracking  phase 
boimdaries  and  intrinsically  includes  all  physical  processes,  material  properties  and  equations  of 
state  in  both  subcritical  and  supercritical  environments.  MD  also  allows  the  tracking  of  injected 
liquid  molecules  into  the  chamber  gaseous  environment.  Simulations  of  subcritical  liquid 
nitrogen  injection  and  atomization  resulted  in  droplet  sizes  that  matched  Rayleigh  theory. 
Subcritical  breakup  length  as  a  function  of  pressure  displayed  behavior  that  matched 
experiments.  Axial  and  radial  density  profiles  for  supercritical  jet  injection  displayed  the  smooth 
transition  of  the  density  profiles  from  the  high  jet  centerline  to  ambient  chamber  conditions  and 
jet  divergence  angles  that  agreed  in  magnitude  with  steady  state  empirical  models. 

Results 


a)  Simulation  of  liquid  nitrogen  injection 

A  molecular  dynamics  simulation  of  liquid  nitrogen  injection  into  gaseous  nitrogen  was 
constructed.  The  nitrogen  molecules  were  simulated  through  the  use  of  constraint  dynamics*’  . 
Each  nitrogen  atom  in  the  molecule  is  simulated  separately  with  its  own  Lennard-Jones  12-6 
interatomic  potential.  A  constraint  force  is  added  to  the  sum  of  the  interatomic  forces  acting  on 
each  atom  and  iteratively  adjusted  using  the  RATTLE  algorithm  to  satisfy  the  specified 
constraint,  which  in  the  case  of  molecular  nihogen  is  the  fixed  bond  length^. 

The  code  is  a  parallel  processor  domain  decomposition  simulation  developed  using  the 
MPI  parallel  library.  The  domain  decomposition  algorithm  divides  the  computational  region 
along  the  longitudinal  axis  into  disk-shaped  geometric  domains,  with  each  domain  assigned  to  a 
separate  processor.  Particles  are  exchanged  between  processors  when  they  travel  between 
domains.  A  domain  decomposition  algorithm  is  the  most  difficult  to  develop  but  is  the  most 
computational  efficient  when  scaled  to  a  large  number  of  particles  and  processors.  A  Verlet 
nearest  neighbor  list^  was  used  to  reduce  the  needed  force  calculations  to  those  atoms  within  a 
cutoff  radius  of  2.5  a . 

Approximately  5000  nitrogen  molecules  were  used  in  the  initial  simulation  region  that 
was  4  nm  in  diameter  and  40  nm  long.  The  atoms  were  initialized  at  a  liquid  density  in  a  cubic 
lattice  with  a  Maxwellian  velocity  distribution  representative  of  a  given  initial  temperature.  A 
stochastic  boundary  condition  was  applied  on  the  sidewalls  using  diffuse  reflections  to  obtain  a 
no-slip  boundary  at  a  constant  temperature.  Periodic  boundary  conditions  were  used  on  the  ends 
of  the  tube,  if  an  atom  exited  from  one  end  of  the  tube  it  re-entered  at  the  other  end.  The  time 
step  used  was  2.5  femtoseconds  (fs).  10,000  time  steps  were  used  to  allow  the  liquid  argon  to 
equilibrate,  then  a  body  force  was  applied  to  the  atoms  to  induce  a  mean  flow  in  the  tube.  The 


*  Sadus,  R,  Molecular  Simulation  of  Fluids^  Elsevier,  Amsterdam,  1999. 

^  Allan,  M.  P.  and  Tildesley,  D.  J.,  Computer  Simulation  of  Liquids^  Oxford  University  Press,  New  York,  1984. 

^  Andersen,  H.  C.,  “Rattle:  a  ‘Velocity’  Version  of  the  Shake  Algorithm  for  Molecular  Dynamics  Calculations,”  */. 
Comp.  Phys.y  Vol.  52, 1983,  pp.  24-34. 


simulation  was  continued  for  another  350,000  time  steps  to  allow  the  fully  developed  pipe  flow 
to  reach  equilibrium. 

Figure  1  shows  the  radial  velocity  profile  of  the  computed  fully  developed  pipe  flow. 
The  analytically  predicted  parabolic  velocity  profile  for  the  Poiseuille  flow  is  obtained.  The 
Reynolds  number  of  the  flow  is  approximately  5.  The  simulation  requires  approximately  twelve 
hours  to  run  on  sixteen  processors  of  our  Beowulf  PC  cluster. 


Fig.  1  Radial  velocity  profile  for  fiilly  developed  liquid  nitrogen  pipe  flow.  Re  -  5. 

After  the  liquid  argon  pipe  flow  reached  equilibrium,  the  initial  simulation  region  was 
mated  to  a  larger  diameter  (20  nm)  cylindrical  region  40  nm  long  containing  quiescent  gaseous 
nitrogen  at  a  specified  pressure  and  temperature.  The  liquid  nitrogen  was  then  allowed  to  flow 
under  the  influence  of  the  body  force  into  the  gaseous  region.  The  periodic  bormdary  conditions 
were  removed  from  the  upstream  end  of  the  liquid  region  and  the  downstream  end  of  the  gaseous 
region.  The  gaseous  region  had  the  same  sidewall  boundary  condition  as  the  liquid  region,  i.e., 
no-slip  and  a  specified  temperature.  An  additional  eight  processors  were  used  for  the  gaseous 
region  utilizing  the  same  domain  decomposition  scheme  as  for  the  liquid  region.  Figure  2  shows 
a  schematic  of  the  initialization  of  the  jet  injection  simulation. 

b)  Implementation  of  molecular  tagging 

In  order  to  follow  the  mixing  of  the  injected  liquid  nitrogen  into  the  gaseous  nitrogen,  a 
means  to  track  the  initially  liquid  and  initially  gas  molecules  was  implemented  in  the  simulation. 
Figure  3  shows  an  initialized  simulation  where  the  molecules  initially  in  the  liquid  are  colored 
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Fig.  2  Schematic  of  initialization  of  molecular  dynamics  jet  injection  simulation. 


Fig.  3  Molecular  tagging:  red  molecules  are  initially  in  the  liquid  while  green  molecules 
are  initially  in  the  gas.  Both  molecules  are  nitrogen. 


red  while  the  molecules  initially  in  the  gas  are  colored  green.  The  molecules  retain  their  initial 
color  as  the  simulation  proceeds,  allowing  the  mixing  of  the  liquid  into  the  gas  to  be  followed. 

c)  Subcritical  liquid  nitrogen  injection 

In  order  the  validate  the  accuracy  of  the  molecular  dynamics  simulation  to  model 
supercritical  injection  and  mixing,  simulations  were  initially  conducted  to  model  subcritical 
liquid  nitrogen  injection  for  comparison  to  classical  analytical  theories.  Figure  4  shows  the 
simulation  of  subcritical  liquid  nitrogen  at  84  K  and  0.5  MPa  injected  into  gaseous  nitrogen  at 
the  same  conditions.  The  blue  molecules  are  initially  liquid  and  the  red  molecules  are  initially 
gas.  The  injection  velocity  was  44  m/s,  corresponding  to  a  Reynolds  number  Re=1.14.  The 
pinch-off  and  formation  of  two  liquid  nitrogen  droplets  can  be  seen.  The  diameter  of  the  drops 
matches  the  Rayleigh  theory  predicted  size  for  low  Reynolds  number  flow  of  1.89  times  the 
diameter  of  the  liquid  jet'*.  Clearly  the  instability  of  the  liquid  jet  due  to  the  surface  tension  is 
reproduced. 


Fig.  4  Molecular  dynamics  simulation  of  subcritical  liquid  nitrogen  injection  at  0.5  MPa  and  84 

K  demonstrating  Rayleigh  break-up. 

Simulations  have  also  been  conducted  at  1, 2  and  3  MPa  liquid  and  gas  pressures  (critical 
pressure  for  nitrogen  is  3.38  MPa).  The  Rayleigh  atomization  is  still  reproduced  at  these  higher 
pressures,  although  the  droplets  are  more  difficult  to  visualize  due  to  the  increased  gas  density. 
The  liquid  core  breakup  length  as  a  fimction  of  pressure  is  shown  in  Figure  5.  The  expected 
decrease  in  breakup  length  with  an  increase  in  ambient  pressure  is  observed^.  The  difference  in 
breakup  length  between  the  higher  ambient  pressme  cases  is  also  less,  indicative  of  higher 
pressures  having  less  of  an  effect  on  breakup  length. 


*  Rayleigh,  Lord,  “On  the  Instability  of  Jets,”  Proc.  Lond.  Math.  Soc.,  Vol.  10, 1878,  pp.  4-13. 

*  Lefebvre,  A.  H.,  Atomization  and  Sprays,  Hemisphere  Publishing,  New  York,  1988,  pp.  57-59. 


Ambient  Gas  Pressure  (MPa) 

Fig.  5  Influence  of  ambient  gas  pressure  on  subcritical  jet  breakup  length  from  molecular 

dynamics  simulations. 

d)  Supercritical  liquid  nitrogen  injection 

Simulations  have  been  initiated  of  liquid  nitrogen  injection  into  a  supercritical  nitrogen 
environment.  The  low  Reynolds  number  of  the  molecular  dynamics  simulations  results  in  a 
laminar  flow  injection  of  the  liquid  nitrogen.  Although  these  can  be  compared  to  measurements 
of  jet  spreading  angle  and  radial  and  axial  density  profiles  taken  from  Raman  measurements  at 
DLR®  and  AFRL’  of  turbulent  liquid  nitrogen  injection,  current  work  at  DLR  is  focussed  on 
repeating  these  measurements  for  laminar  flow  injection  for  comparison  to  the  simulations. 
Previous  work  at  DLR  as  shown  that  the  nondimensionalized  density  profiles  are  not  a  function 
of  the  jet  diameter®,  thus  allowing  the  comparison  of  the  experimental  results  to  the  molecular 
dynamics  simulations.  Figure  6  shows  a  shadowgraph  image  of  a  turbulent  supercritical  nitrogen 
jet  at  4  MPa  and  130  K  taken  at  DLR^.  In  order  to  quantitatively  compare  results  from  the 
molecular  dynamics  simulations  to  the  experimental  measurements,  means  for  calculating  the 
radial  and  axial  density  profiles  produced  by  the  MD  simulations  are  being  developed.  Figure  7 
shows  radial  density  profiles  taken  from  Raman  measurements  for  a  variety  of  pressures  and 
temperatures  from  the  experiments  at  DLR^.  Figure  8  shows  the  radial  density  profile 
progression  at  various  axial  locations  downstream  of  the  injector  obtained  from  a  molecular 
dynamics  simulation  of  nitrogen  injection  at  40  m/s,  4  MPa  and  123  K.  The  smooth  curves  are 
indicative  of  the  supercritical  nature  of  the  flow.  Since  there  is  no  jet  surface,  the  expected 
density  profile  is  a  smooth  transition  from  the  high  jet  centerline  density  dropping  off  to  the 
ambient  chamber  conditions. 


*  Mayer,  W.,  Telaar,  J.,  Branam,  R.,  Schneider,  G.  and  Hussong,  J.,  “Raman  Measurements  of  Cryogenic  Injection 
at  Supercritical  Pressure,”  scheduled  for  publication  in  Int.  J.  of  Heat  and  Mass  Transfer,  2003. 

’  Cheroudi,  B.,  Cohn,  R.  and  Talley,  D.,  “Cryogenic  Shear  Layers:  Experiments  and  Phenomenological  Modeling  of 
the  Initial  Growth  Rate  under  Subcritical  and  Supercritical  Conditions,”  Int.  J.  of  Heat  and  Fluid  Flow,  Vol.  23, 
2002,  pp.  554-563. 

*  Branam,  R.  and  Mayer,  W.,  “Characterization  of  Cryogenic  Injection  at  Supercritical  Pressure,”  Journal  of 
Propulsion  and  Power,  Vol.  19,  No.  3,  May-June  2003,  pp.  342-355. 
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Fig.  6  Shadowgraph  image  of  supercritical  liquid  nitrogen  jet  at  4  MPa  and  130 
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Fig.  7  Raman  measured  nondimensionalized  radial  density  profiles  for  liquid  nitrogen 


injection  . 
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Fig.  8  Radial  density  profile  progression  fi’om  molecular  dynamics  simulation;  4  MPa,  123  K 
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